library(StereoMorph)
library(tidyverse)
Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages ----------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats
library(geomorph)
Loading required package: rgl
Loading required package: ape
library(listviewer)
library(stringr)
(localities<-read_csv("Marks_Mox - Sheet1.csv"))
Parsed with column specification:
cols(
  CatalogNumber = col_integer(),
  ScientificName = col_character(),
  `# of individuals` = col_character(),
  Locality = col_character(),
  Country = col_character(),
  YearCollected = col_integer(),
  Photo = col_character(),
  species = col_character(),
  drainage = col_character(),
  river = col_character()
)
localities %>%
  select(CatalogNumber,species,drainage,river)%>%
  rename(cat_num=CatalogNumber)->loc
#files<-paste0("Mox_images/shapes/TU198404_",1:6,"_L.txt")
a<-readShapes(file = "Mox_images/shapes/",fields=c("landmarks.scaled","curves.scaled"))

Put landmark data into a tibble.

data<-tibble(num=attr(a$landmarks.scaled,"dimnames")[[3]],
       fixedlm=array_branch(a$landmarks.scaled,margin = 3),
       c_body_ant=map(a$curves.scaled,"body_ant"),
       c_body_post=map(a$curves.scaled,"body_post"),
       c_opercle=map(a$curves.scaled,"opercle"))
data %>%
  mutate(cat_num=str_extract(num,"[0-9]+"))->data
data$cat_num<-as.integer(data$cat_num)
(data<-left_join(data,loc))
Joining, by = "cat_num"

convert curves into landmarks, evenly spaced along curves. Bind fixed and semi-landmarks together and remove duplicates.

data %>%
  mutate(body_ant=map(c_body_ant,~pointsAtEvenSpacing(.x,n=10)))%>%
  mutate(body_post=map(c_body_post,~pointsAtEvenSpacing(.x,n=10)))%>%
  mutate(opercle=map(c_opercle,~pointsAtEvenSpacing(.x,n=5))) %>%
  mutate(land_marks=pmap(list(fixedlm,body_ant,body_post,opercle),rbind))%>%
  mutate(land_marks=map(land_marks,~unique(.x)))->data

Convert list (and bind several arrays together) using sapply()

new_a<-sapply(data$land_marks, I, simplify="array")

Generate plot to aid in defining sliding, semi-landmarks using AUTO mode of define.sliders(). There are 20 fixed landmarks, curves are found between landmarks 1 and 2 (anterior, dorsal body), 3 and 4 (posterior, dorsal body), and 10 and 11 (opercle)

dd<-as.data.frame(new_a[,,1])
dd$label<-1:length(dd$V1)
ggplot(dd,aes(V1,V2))+
  geom_point(alpha=0.7)+
  geom_text(label=dd$label,check_overlap = F,nudge_x = 1,size=3)

Generate semi-landmarks matrix for gpagen using define.sliders().

curves<-rbind(define.sliders(c(1,29:36,2)),
              define.sliders(c(3,21:28,4)),
              define.sliders(c(10,37:39,11)))
gpa<-gpagen(new_a,curves = curves)

  |                                                                             
  |                                                                       |   0%
  |                                                                             
  |==============                                                         |  20%
  |                                                                             
  |============================                                           |  40%
  |                                                                             
  |===========================================                            |  60%
  |                                                                             
  |=========================================================              |  80%
  |                                                                             
  |=======================================================================| 100%
gpa

Call:
gpagen(A = new_a, curves = curves) 



Generalized Procrustes Analysis
with Partial Procrustes Superimposition

20 fixed landmarks
19 semilandmarks (sliders)
2-dimensional landmarks
6 GPA iterations to converge
Minimized squared Procrustes Distance used


Consensus (mean) Configuration

                 X             Y
 [1,] -0.189292837 -0.0188987033
 [2,]  0.029581528  0.0433555969
 [3,]  0.103513405  0.0240097538
 [4,]  0.291057191  0.0149602944
 [5,]  0.292731582 -0.0399869752
 [6,]  0.222333149 -0.0435116367
 [7,]  0.181608223 -0.0589950374
 [8,]  0.055101051 -0.0671944886
 [9,] -0.083371192 -0.0516103724
[10,] -0.094910572 -0.0507079052
[11,] -0.101084694  0.0135698686
[12,] -0.147023141 -0.0044754231
[13,] -0.157494327 -0.0047139364
[14,] -0.154682252  0.0023113931
[15,] -0.146676471  0.0066572908
[16,] -0.138713301  0.0041156807
[17,] -0.136003081 -0.0041973511
[18,] -0.138979090 -0.0116328834
[19,] -0.147209135 -0.0150931918
[20,] -0.154832756 -0.0115755409
[21,] -0.172234467  0.0001636571
[22,] -0.150291550  0.0134418372
[23,] -0.126774867  0.0226499323
[24,] -0.102904749  0.0295011512
[25,] -0.077958549  0.0350546513
[26,] -0.051693480  0.0387790044
[27,] -0.024679231  0.0411041028
[28,]  0.003760085  0.0422855421
[29,]  0.124162947  0.0208084767
[30,]  0.145958581  0.0178719752
[31,]  0.167119641  0.0151792999
[32,]  0.188043939  0.0127899438
[33,]  0.208877947  0.0106387953
[34,]  0.229624137  0.0089623537
[35,]  0.250128884  0.0083247394
[36,]  0.270922320  0.0107404774
[37,] -0.087803161 -0.0356431614
[38,] -0.087333514 -0.0180416327
[39,] -0.092578192 -0.0009975782
plotAllSpecimens(gpa$coords,mean=F)

plot PCA

PCA <- plotTangentSpace(gpa$coords,warpgrids = F)

ggplot(pca, aes(PC1, PC2, color = id2)) +
  geom_polygon(data = hull,aes(PC1,PC2,group=id,fill=id),alpha=0.3) +
  geom_point()+
  facet_wrap(~id2)

M<-mshape(gpa$coords)
PC<-PCA$pc.scores[,1:2]
preds<-shape.predictor(gpa$coords,x=PC,Intercept = FALSE,
                       pred1=c(-0.05,0.04)) 
                       
GP<-gridPar(pt.size=0.5,tar.pt.size=0.5,n.col.cell = 50)      
plotRefToTarget(M,preds$pred1,mag = 2,method = "vector",gridPars = GP)

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShTdGVyZW9Nb3JwaCkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2VvbW9ycGgpCmxpYnJhcnkobGlzdHZpZXdlcikKbGlicmFyeShzdHJpbmdyKQpgYGAKCmBgYHtyfQoobG9jYWxpdGllczwtcmVhZF9jc3YoIk1hcmtzX01veCAtIFNoZWV0MS5jc3YiKSkKYGBgCmBgYHtyfQpsb2NhbGl0aWVzICU+JQogIHNlbGVjdChDYXRhbG9nTnVtYmVyLHNwZWNpZXMsZHJhaW5hZ2Uscml2ZXIpJT4lCiAgcmVuYW1lKGNhdF9udW09Q2F0YWxvZ051bWJlciktPmxvYwoKYGBgCgoKYGBge3J9CiNmaWxlczwtcGFzdGUwKCJNb3hfaW1hZ2VzL3NoYXBlcy9UVTE5ODQwNF8iLDE6NiwiX0wudHh0IikKYTwtcmVhZFNoYXBlcyhmaWxlID0gIk1veF9pbWFnZXMvc2hhcGVzLyIsZmllbGRzPWMoImxhbmRtYXJrcy5zY2FsZWQiLCJjdXJ2ZXMuc2NhbGVkIikpCmBgYAoKUHV0IGxhbmRtYXJrIGRhdGEgaW50byBhIHRpYmJsZS4gCmBgYHtyfQpkYXRhPC10aWJibGUobnVtPWF0dHIoYSRsYW5kbWFya3Muc2NhbGVkLCJkaW1uYW1lcyIpW1szXV0sCiAgICAgICBmaXhlZGxtPWFycmF5X2JyYW5jaChhJGxhbmRtYXJrcy5zY2FsZWQsbWFyZ2luID0gMyksCiAgICAgICBjX2JvZHlfYW50PW1hcChhJGN1cnZlcy5zY2FsZWQsImJvZHlfYW50IiksCiAgICAgICBjX2JvZHlfcG9zdD1tYXAoYSRjdXJ2ZXMuc2NhbGVkLCJib2R5X3Bvc3QiKSwKICAgICAgIGNfb3BlcmNsZT1tYXAoYSRjdXJ2ZXMuc2NhbGVkLCJvcGVyY2xlIikpCgpgYGAKCmBgYHtyfQpkYXRhICU+JQogIG11dGF0ZShjYXRfbnVtPXN0cl9leHRyYWN0KG51bSwiWzAtOV0rIikpLT5kYXRhCmBgYAoKYGBge3J9CmRhdGEkY2F0X251bTwtYXMuaW50ZWdlcihkYXRhJGNhdF9udW0pCihkYXRhPC1sZWZ0X2pvaW4oZGF0YSxsb2MpKQpgYGAKCmNvbnZlcnQgY3VydmVzIGludG8gbGFuZG1hcmtzLCBldmVubHkgc3BhY2VkIGFsb25nIGN1cnZlcy4gQmluZCBmaXhlZCBhbmQgc2VtaS1sYW5kbWFya3MgdG9nZXRoZXIgYW5kIHJlbW92ZSBkdXBsaWNhdGVzLiAKYGBge3J9CmRhdGEgJT4lCiAgbXV0YXRlKGJvZHlfYW50PW1hcChjX2JvZHlfYW50LH5wb2ludHNBdEV2ZW5TcGFjaW5nKC54LG49MTApKSklPiUKICBtdXRhdGUoYm9keV9wb3N0PW1hcChjX2JvZHlfcG9zdCx+cG9pbnRzQXRFdmVuU3BhY2luZygueCxuPTEwKSkpJT4lCiAgbXV0YXRlKG9wZXJjbGU9bWFwKGNfb3BlcmNsZSx+cG9pbnRzQXRFdmVuU3BhY2luZygueCxuPTUpKSkgJT4lCiAgbXV0YXRlKGxhbmRfbWFya3M9cG1hcChsaXN0KGZpeGVkbG0sYm9keV9hbnQsYm9keV9wb3N0LG9wZXJjbGUpLHJiaW5kKSklPiUKICBtdXRhdGUobGFuZF9tYXJrcz1tYXAobGFuZF9tYXJrcyx+dW5pcXVlKC54KSkpLT5kYXRhCmBgYAoKQ29udmVydCBsaXN0IChhbmQgYmluZCBzZXZlcmFsIGFycmF5cyB0b2dldGhlcikgdXNpbmcgYHNhcHBseSgpYApgYGB7cn0KbmV3X2E8LXNhcHBseShkYXRhJGxhbmRfbWFya3MsIEksIHNpbXBsaWZ5PSJhcnJheSIpCmBgYAoKCkdlbmVyYXRlIHBsb3QgdG8gYWlkIGluIGRlZmluaW5nIHNsaWRpbmcsIHNlbWktbGFuZG1hcmtzIHVzaW5nIEFVVE8gbW9kZSBvZiBgZGVmaW5lLnNsaWRlcnMoKWAuIFRoZXJlIGFyZSAyMCBmaXhlZCBsYW5kbWFya3MsIGN1cnZlcyBhcmUgZm91bmQgYmV0d2VlbiBsYW5kbWFya3MgMSBhbmQgMiAoYW50ZXJpb3IsIGRvcnNhbCBib2R5KSwgMyBhbmQgNCAocG9zdGVyaW9yLCBkb3JzYWwgYm9keSksIGFuZCAxMCBhbmQgMTEgKG9wZXJjbGUpIAoKYGBge3J9CmRkPC1hcy5kYXRhLmZyYW1lKG5ld19hWywsMV0pCmRkJGxhYmVsPC0xOmxlbmd0aChkZCRWMSkKCmdncGxvdChkZCxhZXMoVjEsVjIpKSsKICBnZW9tX3BvaW50KGFscGhhPTAuNykrCiAgZ2VvbV90ZXh0KGxhYmVsPWRkJGxhYmVsLGNoZWNrX292ZXJsYXAgPSBGLG51ZGdlX3ggPSAxLHNpemU9MykKYGBgCgpHZW5lcmF0ZSBzZW1pLWxhbmRtYXJrcyBtYXRyaXggZm9yIGdwYWdlbiB1c2luZyBgZGVmaW5lLnNsaWRlcnMoKWAuIApgYGB7cn0KY3VydmVzPC1yYmluZChkZWZpbmUuc2xpZGVycyhjKDEsMjk6MzYsMikpLAogICAgICAgICAgICAgIGRlZmluZS5zbGlkZXJzKGMoMywyMToyOCw0KSksCiAgICAgICAgICAgICAgZGVmaW5lLnNsaWRlcnMoYygxMCwzNzozOSwxMSkpKQpgYGAKCmBgYHtyfQpncGE8LWdwYWdlbihuZXdfYSxjdXJ2ZXMgPSBjdXJ2ZXMpCmdwYQpgYGAKYGBge3J9CnBsb3RBbGxTcGVjaW1lbnMoZ3BhJGNvb3JkcyxtZWFuPUYpCmBgYApwbG90IFBDQQpgYGB7cn0KUENBIDwtIHBsb3RUYW5nZW50U3BhY2UoZ3BhJGNvb3Jkcyx3YXJwZ3JpZHMgPSBGKQpgYGAKYGBge3J9CnBjYTwtYXNfdGliYmxlKFBDQSRwYy5zY29yZXNbLGMoMSwyKV0pCnBjYSRpZDwtZGF0YSRkcmFpbmFnZQpwY2EkaWQyPC1kYXRhJHNwZWNpZXMKCgpwY2EgJT4lCiAgc3BsaXQoLiRpZCklPiUKICBtYXAofi54W2NodWxsKHg9LngkUEMxLHk9LngkUEMyKSxdKSAlPiUKICBiaW5kX3Jvd3MoLikgLT4gaHVsbApgYGAKCgoKYGBge3J9CmdncGxvdChwY2EsIGFlcyhQQzEsIFBDMiwgY29sb3IgPSBpZDIpKSArCiAgZ2VvbV9wb2x5Z29uKGRhdGEgPSBodWxsLGFlcyhQQzEsUEMyLGdyb3VwPWlkLGZpbGw9aWQpLGFscGhhPTAuMykgKwogIGdlb21fcG9pbnQoKSsKICBmYWNldF93cmFwKH5pZDIpCmBgYAoKYGBge3J9Ck08LW1zaGFwZShncGEkY29vcmRzKQpQQzwtUENBJHBjLnNjb3Jlc1ssMToyXQpwcmVkczwtc2hhcGUucHJlZGljdG9yKGdwYSRjb29yZHMseD1QQyxJbnRlcmNlcHQgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICBwcmVkMT1jKC0wLjA1LDAuMDQpKSAKICAgICAgICAgICAgICAgICAgICAgICAKR1A8LWdyaWRQYXIocHQuc2l6ZT0wLjUsdGFyLnB0LnNpemU9MC41LG4uY29sLmNlbGwgPSA1MCkgICAgICAKcGxvdFJlZlRvVGFyZ2V0KE0scHJlZHMkcHJlZDEsbWFnID0gMixtZXRob2QgPSAidmVjdG9yIixncmlkUGFycyA9IEdQKQpgYGAKCgoK